A mathematical model of intercellular signaling during epithelial wound healing 
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Recent experiments in epithelial wound healing have demonstrated the necessity of Mitogen- 
activated protein kinase (MAPK) for coordinated cell movement after damage. This MAPK activity 
is characterized by two wave-like phenomena. One MAPK "wave" that originates immediately after 
injury, propagates deep into the cell layer, and then rebounds back to the wound interface. After 
this initial MAPK activity has largely disappeared, a second MAPK front propagates slowly from 
the wound interface and continues into the tissue, maintaining a sustained level of MAPK activity 
throughout the cell layer. It has been suggested that the first wave is initiated by reactive oxygen 
species (ROS) generated at the time of injury. In this paper, we develop a minimal mechanistic 
diffusion-convection model that reproduces the observed behavior. The main ingredients of our 
model are a competition between ligand (e.g., Epithelial Growth Factor) and ROS for the activation 
of Epithelial Growth Factor Receptor (EGFR) and a second MAPK wave that is sustained by stresses 
induced by the slow cell movement that closes the wound. We explore the mathematical properties 
of the model in connection with the bistability of the MAPK cascade and look for traveling wave 
solutions consistent with the experimentally observed MAPK activity patterns. 

PACS numbers: 87.10.Ed, 87.17.Aa. 



Introduction 



Coordinated cell movement is an essential feature of 
many biological processes, such as wound healing, embry- 
onic morphogenesis, and tumor growth^. In wound heal- 
ing, cell migration and cell contraction are the two main 
mechanisms responsible for wound closure. Cell contrac- 
tion is the dominant mechanism in the closing of small 
wounds through the so called "purse-string" process^. 
For larger wounds, cell contraction is not sufficient, and 
surrounding cells must migrate to close larger wounds. 
While the two mechanisms are not mutually exclusive, 
there are cases where cell migration is the only healing 
process^, such as when a strip of cells from an epithelial 
layer is removed^. Despite the existence of experimen- 
tal assays targeting cell migration during wound healing, 
there are still many open questions. For instance, be- 
fore injury, the cells are resting, but after wounding they 
become motile. What mechanical and biochemical phe- 
nomena regulate motility? Is it the availability of free 
space that leads cells to move toward wound closure? 
What determines the speed of cell migration? To be able 
to answer these questions, we need to understand both 
the mechanical and biochemical aspects of cell migra- 
tion, and how they might regulate each other. While the 
physical mechanisms of cell movement have been well- 
studied2£& 7 -, the complex regulation of the wound heal- 
ing process by biochemical signals and feedback pathways 
remains poorly understood. Recent experimental inves- 
tigations of wounds in epithelial tissue have highlighted 
novel properties of the intercellular signaling necessary 
for healing^. 



Matsubayashi et al. (2004) analyzed wounded epithe- 
lial monolayers of Madin-Darby canine kidney (MDCK) 
cells and showed coordinated movement not only by the 
cells at the wound edge, but by several rows of adja- 
cent cells. In their experiments, cell proliferation did 
not play a significant role. Their study also showed that 
the cellular response of the epithelial monolayer is qual- 
itatively characterized by two wave-like activation pat- 
terns of ERK 1/2, a Mitogen Activated Protein Kinase 
(MAPK). These "waves" are characterized by a time- 
dependent front of higher ERK 1/2 concentration that 
initiates at the wound edge and spreads into the cell layer. 
The first rebounding wave-front propagates into the tis- 
sue and back to the wound quickly, the second wave- 
front is slow and sustains MAPK activity in the tissue 
until wound closure. This second, final wave appears to 
be related to cell migration through a positive feedback 
loop, since inhibition of the second MAPK wave halts co- 
ordinated cell movement. Moreover, during this second 
"wave", MAPK is inactive around the healing edges of 
the injured layer, but still active in the migrating cells 
around the open wound*. 

Nikolic et al. (2006) extended the experimental anal- 
ysis in£ by probing the epithelial wound healing assay 
with a novel wound generating technique. They used 
polydimcthylsiloxanc (PDMS) slabs to create two differ- 
ent wounding protocols. A peel-off injury is created by 
growing part of the epithelial monolayer over a PDMS 
slab. Once the slab is peeled-off, taking with it the cells 
grown over the PDMS, it creates a wound that breaks 
some cells at the wound edge, but leaves the space sur- 
rounding the wound free of cellular debris. The other 
protocol consists of the removal of a PDMS slab that 
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forms the boundary of the cellular monolayer, its re- 
moval leaves intact the cells in direct contact with the 
slab, but now provides free space at one edge of the 
monolayer. This "unconstraining" protocol avoids cell 
tearing upon PDMS removal. Experimenting with these 
two techniques, and with the standard scratch wound as- 
say, the authors confirmed the presence of two separate 
waves of MAPK in case of injury. The unconstraining 
experiments expressed only the slow, final wave. In this 
case, cell movement was limited and random, suggesting 
that availability of free space is not enough to generate 
an organized migration of the epithelial layer and that 
mechanical injury is necessary. Accordingly, scratch and 
peel-off experiments resulted in both waves of MAPK 
activation and in collective monolayer migration. Using 
immunofluorescence, Nikolic et al. (2006) were also able 
to identify reactive oxygen species (ROS) as one of the 
key components of the intercellular signaling responsible 
for the observed pattern of MAPK activation. They dis- 
covered that ROS are generated immediately after injury 
and remain present around the wound edge for at least 
the duration of the first rebounding wave. Inhibition of 
ROS by N-acetyl-L-cysteine resulted in the absence of 
both MAPK waves and cell migration. The results from 
the two studies^ are summarized in Table [fl 

In this paper, our aim is to build a mathematical 
model that reproduces an that can be used to analyze 
the properties of MAPK activation during wound heal- 
ing experiments. Because of its relevance and ubiquity, 
the MAPK pathway remains the subject of many com- 
putational and mathematical modeling studies 1 ^. The 
MAPK cascade is a signal transduction pathway that re- 
lays an external stimulus to the cell, it is characterized 
by the sequential activation of three protein kinases 1 ^. 
Huang and Ferrelli^ proposed a system of 18 differen- 
tial equations representing the 10 reactions that com- 
pose the three-kinase MAPK cascade. A mathematical 
and computational analysis of the system showed that 
the cascade has the effect of amplifying an input signal 
(e.g., receptor phosphorylation) in such a way that its 
overall behavior can be compared to that of a coopera- 
tive enzyme^ 3 -. These first computational studies sparked 
additional modeling efforts that, coupled with ongoing 
discoveries by experimentalists, have led to many more 
advanced models that exhibit characteristics (e.g., bista- 
bility, ultrasensitivity, oscillations, etc.) of the MAPK 
signaling pathwa y 14 ' 15 ' 16 ' 17 . 

Here, we are not interested in the intracellular dy- 
namics of MAPK, but rather on its role within the 
wound healing signaling network. For this reason, we 
will treat the MAPK cascade as a black box, using the 
results i n 11 ' 18 to essentially represent the whole cascade 
as a switch for signal transmission. In this approach, 
MAPK/ERK is the output of the "signaling switch". 
ROS are one input that can activate the switch, but 
arc unlikely to be the only one because of the different 
properties of the two activation waves. As suggested in 9 - 
and in other wound healing experiment o 19 ' 20 , other in- 



Experiment 


Results 


Scratch wound 


Two MAPK waves, cell migration 


Closed wound 


No MAPK activity, no cell 
movement 


MAPK inhibition 


No MAPK waves, no cell migration 


Slow wave inhibition 


No cell migration 


Peel-off wound 


Two MAPK waves, cell migration 


Unconstraining 


Slow, second wave only, no cell 
migration 


ROS inhibition 


No MAPK waves, no cell migration 



TABLE I: Summary of the experimental results from«» 



puts are diffusible ligands and their cell receptors. For 
instance, Epidermal Growth Factor (EGF) and EGF Re- 
ceptor (EGFR) play essential roles in promoting cell mi- 
gration, proliferation and wound closur o 21 ' 22 . Moreover, 
positive and negative feedbacks between EGFR signaling 
and the MAPK cascade have been demonstrated exper- 
imentally and verified computationall y 23 ' 24 ' 25 . Although 
Nikolic et al. suggest EGF as a possible signal, as well 
as other molecules, they did not pursue the topic in 
their work. However, they did identify Reactive Oxy- 
gen Species (ROS) as direct regulator of MAPK activity 
in their wound healing experiments. Since ROS have 
been shown to induce EGFR activation in the absence 
of EGF—, this finding is in agreement with other stud- 
ies that showed the regulatory role of ROS in wound 
healing 2 !^ and MAPK signalingS 9 .^. 

Because of the documented connection between dif- 
fusible signals, wound healing and MAPK activation, 
we propose a mechanistic model based on ligand- 
mediated intercellular signaling that reproduces the ob- 
served MAPK activation pattern, and that is consistent 
with the qualitative experimental features listed in Table 



Mathematical Model 



The experimental results from 8 ^ 9 - provide evidence of 
spatio-temporal MAPK signaling for the regulation of 
cell migration during wound healing, without determin- 
ing the exact biochemical events that govern it. The 
fact that free space by itself is not enough to produce 
coordinated cell migration suggests that ROS is not the 
only diffusible signal needed to generate the two waves 
of MAPK activation. To explain the observed spatio- 
temporal MAPK pattern, at least two diffusible signals 
are needed. Both ROS and EGF molecules are able to dif- 
fuse in the extracellular space (ROS can also move across 
the cell membrane) and both can phosphorylate the EGF 
membrane receptor (EGFR), activating the MAPK cas- 
cade. EGF induces phosphorylation of the cytoplasmic 
tail of EGFR by binding to it. Reynolds et al. (2003) 
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FIG. 1: Schematic of biophysical events during wound heal- 
ing, (a) diffusible ligands (L) phosphorylate membrane recep- 
tors by binding to them, activating the MAPK cascade that 
leads to the production of intracellular protease (P). Protease 
induces ligand release, (b) ROS (S) can also interact with 
membrane receptors by inducing phosphorylation of their cy- 
toplasmic tail, which also activates the MAPK cascade and 
promotes release of ligands into the extracellular matrix. Dif- 
fusible ligand and ROS represent two independent triggers of 
the MAPK cascade, (c) Schematic of the three-kinase cas- 
cade that characterize MAPK signaling, (d) The stresses 
caused by cell movement toward wound closure can lead to 



showed that ROS can induce EGFR phosporylation even 
in the absence of EGF by binding to intracellular phos- 
phatases. Also well documented is the positive feed- 
back between EGF and the MAPK cascade, and its abil- 
ity to produce long range signaling through autocrine 
relays^. Finally, ROS can be generated by mechanical 



stresses like the ones generated by migration of the ep- 
ithelial monolayer—, thus providing a feedback loop be- 
tween MAPK activation, cell motility, and further ROS 
production. These four signaling mechanisms arc sum- 
marized in Fig. [1] and constitute the foundation of our 
mathematical model. 

The wound healing system is a three-dimensional one, 
and the migration of individual cells toward wound clo- 
sure is not exactly normal to the wound edge as shown by 
cell tracking experiments^. Here, we simplify the anal- 
ysis by considering the cell layer in cross section as a 
semi-infinite straight line, with the wound initially posi- 
tioned at the origin. The resulting two-dimensional sys- 
tem consists of a semi-infinite cell layer that is immersed 
in medium of infinite height (the medium in^ is 3mm 
which is much larger than any other length scale in the 
problem). We model four species: ROS, one diffusible 
ligand {e.g., EGF), one ligand receptor {e.g., EGFR) 
and playing the role of the output of the MAPK cascade 
"black box", a protease that is the intracellular precur- 
sor of the ligand {e.g., the piece completing the feedback 
loop between EGF and the MAPK cascade in Figs, 
(a) and[T]-(b)). We denote the local concentrations of 
these species by L, S, R, and P respectively. As de- 
picted in Fig. [TJ the signal can be transmitted to the cell 
in two different ways: through a ligand-receptor complex 
(C L =R-L) and a ROS-receptor complex^ (C s = R-S). 
We will assume that the number of available cell mem- 
brane receptors is in excess, implying that R is approxi- 
mately constant and that ROS-ligand-rcccptor complexes 
are negligible. 

A schematic representation of the system is given in 
Fig. [21 The governing equations and boundary condi- 
tions of our model are 



dL(X,Z,T) f d*L(X,Z,T) , d 2 L(X,Z,T) \ 

df = ° L { dX* + W* )MX,Z = co t T) = 0, 



Dl l gz = k onRL{X, 0, T) - k^C L (X, T) - g h P(X, T), (2) 

dCh QT T) = £ a RL(X, 0, T) - (fc D L ff + fc c L c ) C h (X, T), (3) 
dS(X,Z,T) f d*S(X,Z,T) 8 2 S{X,Z,T) \ 

df— = Ds { 8X* + 80 )^{X,Z = o C ,T) = 0, (4) 

D s dS(X d f T) = k s on RS{X,0,T)-g s n s {C L ,X,T), (5) 

dCS gf T) = k s on RS{X, 0, T) k s cc C s {X, T), (6) 

dP f T T) = -k P P{X,T)+g P U P {C L ,C s ). (7) 

I 



Equations JT]) and (|U) describe the diffusion of ligands and ROS in the extracellular medium with homogeneous 
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diffusion constant Dl and D$, respectiveiy. Eqn. ([2]) 
accounts for the flux of ligand across the surface of the 
cellular layer including ligand-rcceptor complex forma- 
tion with rate constant k^ n , ligand-receptor complex dis- 
sociation with rate constant k^ s , and extracellular ligand 
release by intracellular protease with rate <?l- Eqn. ([3|) 
governs the kinetics of ligand-receptor complexes, with 
new complexes forming at rate k^ n and dissociating with 
rate k^ s , k^ c represents the rate of receptor- mediated cn- 
docytosis of the ligand-bound receptor complexes. Eqn. 
([5]) describes the diffusive flux of ROS due to formation of 
ROS-receptor complexes with rate k£ n and to ROS pro- 
duction by the functional lis (Cl ,X,T). The kinetics of 
ROS-receptor complexes is represented by Eqn. ((6]) and 
its terms are analogous to the ones in Eqn. ([3]), except 
that there is no ROS release from the ROS-receptor com- 
plex in accordance with^S. The last equation describes 
the cellular response to extracellular signaling through 
the activity of intracellular proteases. Within the wound 
healing framework, protease activity is directly related 
to ERK1/2 activity measured in£. In particular, the pro- 
tease dynamics is characterized by a degradation term 
with rate constant kp and a source term IIp(Cl, Cs) with 
maximum production rate gp. To complete the descrip- 
tion of the mathematical model, we need to impose rea- 
sonable functional forms for ITp and Ilg. 

The role of the functional I1p(Cl,Cs) is to represent 
the intermediate biochemical steps that lead to protease 
production. These steps include the MAPK cascade and 
any other reaction in the feedback loop between ligand 
binding and ligand release {e.g., the solid box in Fig. Q} 
(c)). In the literature, lip is usually represented as a 
sigmoidal function of cell surface complexes such as the 
Hill functio n 11 ! 18 ! 34 . If the level of receptor signaling is 
given by the total concentration of complexes (Cl + Cs), 
we propose the following functional form: 



n P (c L ,c s ) 



(Cl + CsY 



Cl + (C L + C s 



(8) 



where n is an effective Hill coefficient, and Ca. represents 
an activation threshold of the signaling pathway. 

Defining the ROS source ns is more problematic since 
the experimental evidence suggests an interplay between 
cellular signaling and cell migration, thus involving me- 
chanical forces whose description goes beyond the scope 
of this paper. Conversely, the production of ROS due 
to wound induction is embedded in the initial conditions 
of the system and it is not described in ns- Generally, 
ROS production increases with cell metabolism 2 ^ In our 
case, metabolic increase may be related to cells becom- 
ing motil o 29 ' 33 and/or to ligand signaling 3 ^. We use a Hill 
function multiplied by a decaying exponential to repre- 
sent ns: 



U S (C L ,X,T) = 





(C L (X,T-T D )) m 
S%+(Cl(X,T-T d ))-" 



T<T D 
T>T D 



(9) 



The Hill functional represents ligand-mediated ROS pro- 
duction stemming from the phosphorylation of a recep- 
tor's tail during ligand binding. The delay Tp> represents 
the delay between ligand binding and ROS production. 
Finally, the exponential factor in Eqn. ([9]) describes re- 
duction in ROS production due to the decrease in motil- 
ity from cells near the wound edge (X = 0) to cells farther 
from it. 

Upon introducing the following dimcnsionlcss quanti- 
ties 



t = kpT, x = X\ 



I = L- 



kpk^ c k^ n R 

'SLSP ( k o« + k ec. 



z = Z A 



CL = Cl- 



.9L.gp 



cs = C s - 



ft cc 

.9s 



5s 



D kp 

gp 



we express the system of equations in dimensionless form 
dl(x,z,t) d 2 l(x,z,t) d 2 l(x,z,t) 



at 

a l z (x,0,i) 

9c L (a;, t) 
dc s (x,t) 
ds(x, z, t) 

at 



dx 2 



+ 



dz 2 



(10) 



l(x,0,t) - 

(3[l(x,0,t)-CL(x,t)]-p(x,t), (11) 



= V 



l(x,0,t) - c h (x,t), 

s(x, 0, t) — cs(x, t), 

' d 2 s(x, z,t) d 2 s(x,z,t) 



dx 2 



dz< 



v s z (x, 0, t) = s(x,0, t) -tt s (cl,x, t), 



dp 
di 



— = -p + T7 p {c L ,C S ) , 



(12) 

(13) 

(14) 
(15) 
(16) 



where the dimensionless parameters are 



' r h.L ' " 



kp 



V.L I L.L ' 
K oS + K cc 



V 



Dl 




(17) 



The parameters a and v characterize the relative rates of 
diffusion and binding, while P represents the strength of 
complex degradation relative to ligand dissociation. The 
parameters e and S describe the speed of binding and 
endocytosis relative to ligand release mediated by intra- 
cellular species, and r\ is the ratio of diffusivity between 
ROS and EGF ligand. 
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s = 0.1,c A =0.7 



[nucleus ' 



• - S 

• - L 
■ - P 



C 



H 0.5 




FIG. 2: Schematic of simplified signaling pathway during 
wound healing: ligands (L) and ROS (S) are free to dif- 
fuse and bind to cell surface receptors. This binding forms 
complexes (Cs and Cl) that activate MAPK signaling. Com- 
plexes are lost due to endocytosis or, in the case of ligand- 
receptor complexes, through ligand unbinding and endocyto- 
sis. Intracellular proteases (P) release extracellular ligands 
(L). 



The dimensionless protease and ROS production func- 
tions become 



7t p (cl,cs) 



TT s (c L ,X,t) 



(cl + 7 cs)" 



(18) 







t < T 



(c L (x,t-r)) m _ e . Xx t>r ,(19) 



(c L (x,t- t)) 



where c A = C A {kpk^ c )/(g L g P ), 7 = (gskpk^ c )/(g L g P k%), 
t = T D fc P , A = k x L, and s A = (S A kpk^ c )/(g L gp). 



FIG. 3: Steady state configurations, (a) Graphical solution 
of Eqn. (TJHI for s = 0.1, and c A = 0.7. There are three 
solutions representing two stable steady states (Z u and h) and 
one unstable equilibrium (h). (b) If we fix the activation 
threshold (ca = 0.7), different values of ROS concentration 
lead to different steady state configurations. Bistability is 
possible if there is enough ROS, otherwise the only stable 
steady state is the one with no MAPK activity, (c) We fix 
ROS concentration to s = 0.1 and graphically solve Eqn. (|26[) 
for different values of c\. Bi-stability can arise only if the 
activation threshold is sufficiently small. The Hill coefficient 
n — 8 was used in all plots. 



where all functions now represent "outer" solutions valid 
at times beyond initial transients in complex formation. 
We verified that this approximation holds throughout all 
of the analysis performed in the next section. 



Fast Binding Approximation 



Before we proceed to the analysis of the model, we are 
going to make an assumption that significantly simpli- 
fies the model and that is also justifiable biophysically. 
We assume that ligand dissociation and complex degra- 
dation is fast compared to protease degradation, e.g., 
kp <C fc| c , k^ s , k^ c . In this limit, e and S are small and 
Eqns. (fP2")l and (|13[) can be treated as a singular per- 
turbation. On timescales of protease degradation the 
concentration of ligand-receptor and ROS-receptor com- 
plexes are approximately that of the surface concentra- 
tion of free ligand and ROS, respectively. If we consider 
only the "outer" solutions of Eqns. (fH?)) and (flU)) , our 
full model reduces to the three equations: 



Analysis & Results 



The spatio-temporal MAPK activation patterns can 
arise from different mechanisms. For example, one (or 
more) activation pattern could consist of a traveling front 
connecting two stable steady states, corresponding to 
high and low MAPK concentrations. In this section we 
describe the steady states of the system of Eqns. ([20]) - 
(|24p and present an overview of the qualitative behavior 
of the solutions of the model. After establishing the dy- 
namics of the mathematical model, we use the known 
model parameters to determine the nature of the MAPK 
patterns and some of their properties. 



dl(x,z,t) d 2 l(x 1 z 1 t) d 2 l(x,z,t) 



dt 

a l z {x,0,t) 
ds(x, z, t) 
dt 

v s z (x,0,t) 
dp 
dt 



dx 2 



Or. 



l(x,0,t) ~p(x,t), 

'd 2 s(x,z,t) d 2 s(x,z,t) 



n 



dx 2 dz 2 
s(x,0,t) - n a (l,x,t), 

-p + ir p (1,S) , 



(20) 
(21) 

(22) 
(23) 
(24) 



The complexity of our model requires analysis through 
numerical simulations. For this purpose we use an ex- 
plicit finite difference scheme that is forward in time and 
centered in space, implemented in Fortran. We use a 
uniform grid discretization along the direction of the cell 
layer (e.g., x) and a geometrical grid discretization along 
the direction normal to the cell layer (e.g., z) to max- 
imize accuracy and minimize run-time. This approach 
and its advantages have been previously described ir*2^. 
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0.4 H 1 1 1 1 1 1 

0.0 0.1 0.2 0.3 0.4 0.5 0.6 

J 

FIG. 4: Traveling waves connecting the state of MAPK ac- 
tivation and the state of MAPK inactivity can move toward 
or away from the wound depending on ROS concentration at 
cell layer level (s) and activation threshold (ca) when n — > oo. 
The thick black line indicates conditions under which the front 
is not moving. Regimes above the line lead to traveling waves 
moving toward the wound, while those below the line lead to 
fronts moving away from the wound. 

Steady States & Traveling Fronts 

The steady states of the model in Eqs. l|2"0] ) -([2"3 ]) sat- 
isfy: 

I — p, s = 7r s (Z, x, t), and p = ir p (l, s), (25) 

where the overbar indicates that the value of ligand or 
ROS concentration is taken at z = (e.g., I = l(x,z = 
0,t)). The resulting condition 

1=tt p (1,s) (26) 

is always satisfied by the trivial solution (/ = 0, s = 0), 
but under certain conditions it can have two more solu- 
tions as highlighted in Fig. [3]- (a). In this case the three 
roots are two stable steady states, Iq and I2, and an un- 
stable one, l\. The two stable steady states represent a 
state of no ligand signaling, Iq = 0, and a state of active 
ligand signaling, 1% > 0, respectively. From Fig. [3J we 
can also infer how ROS concentration s and the activa- 
tion threshold ca control the steady states of the system. 
If we fix the activation threshold at a sufficiently high 
value, the system attains only the trivial steady state, 
I = Iq = 0, unless there is enough ROS to sustain the 
signaling pathway, as shown in Fig. EJ-(b). Conversely, if 
we fix s, the system will be bistable only if the activation 
threshold ca is sufficiently small (Fig. [3J-(c)). 

Bistability implies that the model admits traveling 
front solutions connecting the two stable steady states. 



We can approximate the front speed analytically by con- 
sidering a simpler scenario in which the concentration 
of ROS at cell layer level, s, is constant and take the 
limit n — ► 00 for the Hill coefficient of 7r p . In this limit, 
the sigmoidal protease production function 7r p equals the 
Hcavisidc function centered at ca — s: 



lim 7r p (Z, a) = H(l + a - c A ) = { i I ~ CA * ■ (27) 

n~*oo I < CA — S 

In this limit the system is bistable for < ca — s < 1, 
and the roots of Eqn. [26] are Iq = 0, l\ = ca — s, and 
I2 = 1. Furthermore, we can determine the velocity and 
direction of the traveling fronts by proceeding as in££, 
obtaining: 



- 1 f ayjq 2 -vq 

TJo, 9(1 + vq)(a z q z - a z vq + 1) 

which gives an implicit relation for v, the velocity of the 
traveling wave for a fixed concentration of ROS. The in- 
tegration variable q in the above integral arises from the 
Fourier transformation used to derive Eqn. 1281 From 
Eqn. [28] we find that the front velocity is a monoton- 
ically increasing function of the parameter a (see Eqn. 
[T7|) . This result is expected since an increase in a cor- 
responds to either an increase in ligand diffusivity or a 
decrease in ligand binding, and both changes result in 
the front reaching farther distances in a shorter amount 
of time. The direction of the front is determined by the 
threshold ca — s. If ca — s < 1/2, the front of active 
MAPK will move away from the wound, and deep into 
the cell layer. If ca — s > 1/2, MAPK activity will recede 
toward the wound. Regimes that delineate forward and 
backward traveling MAPK waves are indicated in Fig. [4] 
These results provide useful insight for the general case of 
n < 00 and diffusing ROS. Numerical simulations show 
that for Hill coefficient as low as n = 6, the instantaneous 
front velocity is within 10% of that obtained from Eqn. 
1281 To summarize, we showed that the system can have 
two stable steady states and that traveling wave solutions 
connecting them are possible. In particular, ROS can de- 
termine the existence, velocity and direction of the fronts 
by effectively regulating the activation threshold of the 
MAPK cascade. 



ROS/EGF regulation of MAPK activation 

Bistability is necessary but not sufficient for the ex- 
istence of traveling wave solutions. In this section we 
explore the parameter space of the wound healing assay 
to determine the nature of the MAPK fronts observed 
To avoid ambiguity, we divide the MAPK dynam- 
ics during wound healing into three wave- like events. The 
first event corresponds to the fast activation of MAPK 
initiated by the wound. It lasts until the activation front 
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Parameter Typical Value Ref. 





10" 8 - 10" 6 cmV 1 


31 


"'on 


10~ 15 - 10~ 12 cm 3 s"' 


31 


k h 

K oB 


1(T 3 - 10 -2 s" 1 


31 


A- L 


10" 3 - 10" 2 s _1 


31 


fcp 


10~ 4 - KT 3 s^ 1 


31 


R 


10 1() - 10 13 cm^ 2 


31 


ffp 


0.17 x 10* cm" 2 s" 1 


30 




0.54 x 10~ 2 s^ 1 


37 


Ca 


10 9 cm -2 


37 



TABLE II: Typical values of model parameters. 



reaches its maximum depth in the cell layer. The sec- 
ond event is characterized by decrease of MAPK activ- 
ity. It moves from deep into the epithelial monolayer 
toward the wound. These first two events qualitatively 
correspond to the first "rebounding wave" observed in 
the experiments^. The last event consists of a slow ac- 
tivation front that starts at the wound edge and moves 
away from the wound. This last "wave" is initiated when 
the cells in the layer start moving to close the wound it- 
self, and is sustained when the wound is large, preventing 
closure in finite time^. 

To reproduce the observed MAPK dynamics, we nu- 
merically integrated Eqns. (|20p - (|24|) using the ligand re- 
lated parameters given in Table HU Although we could 
not find analogous references for physical parameters 
of ROS, we estimated parameter values from various 
sources. We used the self-diffusivity of water together 
with the Einstein relation to bound the value of r) be- 
tween 10 and 100. From the results ir*2£ it seems rea- 
sonable to assume fcf c « k^ c . We assume that the ini- 
tial concentration of ligand and protease is zero, while 
the concentration of ROS is represented by a narrow 
Gaussian with width equal to the size of a single cell; 
it represents the ROS initially released by cell rupture. 
Our numerical results are summarized in Figs. O and [5] 
Fig. [5] compares the time evolution of the distance of the 
front from the wound as predicted by Ecins. (f20]) - (|24p with 
the experimental values observed in^ for scratch wounds. 
The position of the front is determined by evaluating 
the inflection point of protease concentration after each 
time step. The model is able to replicate the observed 
MAPK behavior and we used it to investigate the dy- 
namics of the three activation wave-like patterns. Fig. 
|5] shows the time evolution of the profiles of the three 
MAPK "waves". The first wave (Fig. [6]-(a)) is driven by 
ROS production at the onset of wound and its fast diffu- 
sion. However, there is not enough ROS for either ligand 
or protease concentration to reach the signaling steady 
state (e.g., p = I = 1). As a result, the first "wave" 
can only propagate as far as ~ 480/im before receding. 
As ROS diffuses away, protease concentrations decrease 
(Fig. 0(b)) until the cells start to move (after about 30 



600 n 




min 



FIG. 5: Time evolution of the distance from the wound edge 
of the MAPK front. The circles and the intersecting vertical 
bars represent the average front position and one standard 
deviation error bars, respectively. Their values have been 
obtained from the experimental data irA The solid trace has 
been obtained through a computer simulation of Eqns. (|20[1 - 
(|24[) with dimensionless parameters: a = 0.35, r\ = 16, v — 5, 
c A = 0.625, n = 6, 7 = 2, p A = 0.975, m = 9, A = 0.55, 
t = 2.56. These parameters are based on fcp = 5.6 x 10 _3 s _1 , 
D L = 5.6 x 10- 7 cm 2 s-\ /& = f£ g = 10 _2 s~ 1 , g L = 0.54 x 
10 _2 s- 1 , g P = 0.17 x 10 8 cm- 2 s-\ R = 3.2 x 10 n cm- 2 , k^ n = 
lO^cmV 1 , Ca = 10 9 cm~ 2 . 



minutes from injury). At that time, ROS is produced 
by the moving cells and fuels the positive feedback loop 
between ligand and protease. The nonlinear effects of ir p 
allow protease levels to increase until they reach a sig- 
naling steady-state. At this point, the front moves like 
a true traveling- wave (Fig. El-(c)), with its speed and 
distance traveled regulated by the ROS source function 
7r s . We can also track the time-evolution of the variables 
in the model. Fig. [6]-(d) shows how protease concentra- 
tion at the wound edge x = changes in time. From 
this graph we observe that during the first MAPK event 
the protease concentration never reaches the "signaling" 
steady-state p « 1 and eventually decreases. During the 
third wave, protease concentration reaches the "signal- 
ing" steady-state and remains there as shown by the flat 
part of the graph in Fig. [6]-(d). To summarize, only 
the third MAPK front behaves as a true traveling wave, 
while the initial two events (corresponding to the first 
"rebounding wave" observed in experiments) are actu- 
ally transient, diffusion-driven patterns. 



Conclusions 

We have formulated a mathematical model for the dy- 
namics of intercellular signaling observed during wound 
healing experiments. From this model, we were able 
to replicate the signaling patterns observed in^, and 
to provide insight regarding their nature. Our choice 
of EGF as signaling ligand is based upon literature re- 
view, but lacks experimental evidence within the epithe- 
lial wound healing assay. However, we showed that the 
properties of EGF (and EGFR) fit the profile for the 
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FIG. 6: Time evolution of protease, (a) After injury, the 
fast diffusion of ROS drives MAPK activation, (b) As ROS 
diffuses away, the positive feedback loop between MAPK and 
ligand is not strong enough to sustain signaling and the depth 
of the front decreases, (c) Once the cells start moving, ROS 
production and the ligand-protease feedback loop fuel the ac- 
tivation front, (d) Time evolution of protease concentration 
at x = 0. The parameters used to obtain these plots are the 
same as the ones used to generate Fig. [5] 

unidentified diffusible signals mentioned irA Our model 
can be expanded to incorporate other diffusible signaling 
molecules. Although it may be possible to find physio- 
logically realistic sets of parameters that lead to signaling 
patterns consisting of three separate traveling waves, the 
parameters associated with the EGF/ROS/MAPK sys- 
tem lead to only one final traveling wave. The first two 



notable events being described by purely diffusive and 
decays dynamics, respectively. 

An aspect of our current model that needs improve- 
ment and further analysis is the determination of the 
ROS source function LTg. From the experiments in-, 
this term seems to be negligible since they only detected 
the presence of extracellular ROS up to ~ lOmin 
after wounding. However, there is substantial evidence 
indicating that cell motility and ligand-receptor binding 
can induce ROS production. A plausible explanation 
for these conflicting results could be that ROS produced 
after wounding is fully recaptured by intracellular 
processes (including EGFR phosphorylation) and never 
crosses the cell membrane. Nonetheless, we performed 
many numerical tests and found that if lis = (data not 
shown), then all three MAPK events are diffusion driven 
and the signaling pattern is due to the different diffusion 
properties of ROS and ligand (EGF). A more physically 
realistic approach might be to include the mechanical 
events that take place within the cell layer and the 
reaction that lead to ROS production after ligand bind- 
ing. Such an approach could provide important insights 
about the mechanisms of post -wound ROS production 
and their relevance within the MAPK signaling context. 
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